##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--ccf_file"), type = "character") ,
    make_option(c("--gene_list"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    ccf_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutationTime/result/All_CCF_mutTime.addShare.tsv"
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    gene_list <- "~/20220915_gastric_multiple/dna_combinePublic/images/selectGCClone/GCClone_gene.all_record.list"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/DriverChoose/GeneLOH"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene_list <- opt$gene_list
sample_info <- opt$sample_info
out_path <- opt$out_path
ccf_file <- opt$ccf_file
type <- opt$type

###########################################################################################

dir.create(out_path , recursive = T)
col <- c( "#006699","#DDA520"  )

###########################################################################################

dat_info_all <- data.frame(fread(sample_info))
dat_ccf_all <- data.frame(fread(ccf_file))

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")
geneN <- "TP53"

###########################################################################################
result_use <- c()
for( type in c( "All" , "IGC" , "DGC") ){
    ## 关注的病理亚型
    if(type == "IGC"){
        dat_info <- subset( dat_info_all , Type == "IM + IGC")
    }else if(type == "DGC"){
        dat_info <- subset( dat_info_all , Type == "IM + DGC")
    }else if(type == "All"){
        dat_info <- subset( dat_info_all , Type != "IM + IGC + DGC")
    }

    dat_ccf <- subset( dat_ccf_all , ID %in% dat_info$ID )

    ###########################################################################################
    result_plot <- c()

    tmp <- subset( dat_ccf , Variant_Classification %in% Variant_Type & Hugo_Symbol == geneN )

    if(nrow(tmp) > 0 ){
        tmp$LOH <- ifelse( tmp$minor_cn == 0 & tmp$total_cn==2 , "LOH" , "Non-LOH" )

        tmp <- subset( tmp , Class != "IM" )

        res_tmp <- c()
        for( sample in unique(tmp$ID) ) {
            tmp_use <- subset( tmp , ID == sample )
            tmp_use <- unique(tmp_use[,c("ID" , "CLS" , "LOH")])

            ## 一个人若发生多个突变，算最早的
            if( length(which(tmp_use$CLS=="clonal [share]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [share]" )
            }else if( length(which(tmp_use$CLS=="clonal [early]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [early]" )
            }else if( length(which(tmp_use$CLS=="clonal [NA]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [NA]" )
            }else if( length(which(tmp_use$CLS=="clonal [late]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [late]" )
            }else if( length(which(tmp_use$CLS=="subclonal")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "subclonal" )
            }
            res_tmp <- rbind(res_tmp , tmp_use)
        }

        res_tmp <- data.frame(table(res_tmp$CLS , res_tmp$LOH))

        ## 其余的时间的合并定义为相对较晚
        other_loh <- sum(res_tmp[res_tmp$Var2=="LOH" & !(res_tmp$Var1 %in% c("clonal [share]")),"Freq"])
        other_loh <- data.frame( Var1 = "Non-Share" , Var2 = "LOH" , Freq = other_loh )
        other_nonloh <- sum(res_tmp[res_tmp$Var2=="Non-LOH" & !(res_tmp$Var1 %in% c("clonal [share]")),"Freq"])
        other_nonloh <- data.frame( Var1 = "Non-Share" , Var2 = "Non-LOH" , Freq = other_nonloh )
        other <- rbind( other_loh , other_nonloh ) 

        res_tmp <- rbind( res_tmp , other )
        res_tmp <- res_tmp %>%
        group_by( Var1 ) %>%
        summarize( Var2 = Var2 , Freq = Freq ,  Ratio = Freq/sum(Freq) )
        res_tmp$Hugo_Symbol <- geneN
        result_plot <- rbind(result_plot , res_tmp)
    }

    ###########################################################################################
    ## 基因分布堆叠图
    dat <- result_plot
    dat$Ratio[is.na(dat$Ratio)] <- 0
    dat$value_text <- paste0( round(dat$Ratio , 2) * 100 , "%")
    dat <- subset( dat , Var1 %in%  c("Non-Share" , "clonal [share]") )
    dat$Var1 <- ifelse( dat$Var1 == "clonal [share]" , "Share" , "Non-Share" )

    ## 计算p值
    trans <- function(num){
        up <- floor(log10(num))
        down <- round(num*10^(-up),2)
        text <- paste("p == ",down," %*% 10","^",up)
        return(text)
    }

    result <- c()
    for( geneN in unique(dat$Hugo_Symbol) ){
        tmp <- subset( dat , Hugo_Symbol == geneN )

        a <- subset( tmp , Var1 == "Non-Share" & Var2 == "LOH" )$Freq
        b <- subset( tmp , Var1 == "Non-Share" & Var2 == "Non-LOH" )$Freq
        c <- subset( tmp , Var1 == "Share" & Var2 == "LOH" )$Freq
        d <- subset( tmp , Var1 == "Share" & Var2 == "Non-LOH" )$Freq

        if(length(a)==0){a=0}
        if(length(b)==0){b=0}
        if(length(c)==0){c=0}
        if(length(d)==0){d=0}

        p <- fisher.test( matrix( c(a,b,c,d) , ncol = 2 ) )$p.value

        if( p < 0.001 ){
            p_text <- trans(p)
        }else{
            p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
        }

        tmp$p_text <- ""
        tmp$p_text[1] <- p_text
        result <- rbind( result , tmp )
    }

    result$Var2 <- factor( result$Var2 , levels = c("LOH" , "Non-LOH") )
    result$Var1 <- factor( result$Var1 , levels = c("Share" , "Non-Share") )

    show_gene <- unique(subset( result , Var1=="Share" )$Hugo_Symbol)
    show_gene <- c("TP53" , "APC" , "PIK3CA" , "CDH1")
    show_gene <- c("TP53")
    result <- subset( result , Hugo_Symbol %in% show_gene  )
    #gene_order <- c(
    #    "TP53" , "ARID1A" , "CDH1" , "APC" , 
    #    "ERBB2" , "PIK3CA" , "RNF43" , "MAP2K7" ,
    #    "MTRR" , "MUC6" , "CFTR" , "BMP6" , "GAL3ST3")
    result$Hugo_Symbol <- factor(result$Hugo_Symbol , levels = show_gene )
    result$Type <- type

    result_use <- rbind(result_use , result)
}

plot <- ggplot( data = result_use , aes( x = Var1 , y = Ratio , fill = Var2 ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion")+
    facet_grid(.~Type) +
    theme(panel.grid = element_blank())+
    scale_fill_npg() +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 8,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold',angle = 45, vjust = 1, hjust=1) ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/Driver_Trunk.TP53.LOH.pdf")
ggsave( out_name , plot , width = 6 , height = 4 )
